Adapted and edited by Joseph White

This notebooks will help me talk through a bit of material regarding spatial data, with some introduction to spatial joins.

Key topics for today:

Packages

Standards:

library(knitr)
library(tidyverse)
library(janitor)
library(lubridate) # because we will probably see some dates
library(here) # a package I haven't taught you about before that doesn't do much, but ....
library(rnaturalearth)
library(WDI)
library(tigris)
library(rgdal)
library(sp)

Some additional packages focused on today’s work:

library(sf) # working with simple features - geospatial
library(tmap)
library(tidycensus)

Informational resources

Using the Neighborhood Geospatial Data (using /data)

Our first data source comes from opendata.dc

https://opendata.dc.gov/datasets/DCGIS::dc-health-planning-neighborhoods/about

I will use the GeoJSON file. (Newer, not necessarily better, but … a single file. Not smaller, but … this one is not big.)

Data is easily readable

neigh=st_read(here("DC_Health_Planning_Neighborhoods_joey.geojson")) %>% clean_names()
Reading layer `DC_Health_Planning_Neighborhoods' from data source 
  `/Users/joewhite/Desktop/R-Studio/ds241/ds241_portfolio/DC_Health_Planning_Neighborhoods_joey.geojson' 
  using driver `GeoJSON'
Simple feature collection with 51 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.11976 ymin: 38.79165 xmax: -76.9094 ymax: 38.99556
Geodetic CRS:  WGS 84
class(neigh)
[1] "sf"         "data.frame"
plot(neigh)

Reminder - Joins

df1=tibble(fruit=c("apple","banana","cherry"),cost=c(1.5,1.2,2.25))
df2=tibble(fruit=c("apple","apple","cherry","lemon"),
           desert=c("pie","cobbler","cobbler","cheesecake"),
           cal=c(400,430,500,550))
df1
df2
left_join(df1,df2,by="fruit")

Investigating joining spatial and non-spatial data

Covid case information is available from opendatadc:

https://opendata.dc.gov/datasets/DCGIS::dc-covid-19-total-positive-cases-by-neighborhood/about

Read cases information:

df_c=read_csv(here("DC_COVID-19_Total_Positive_Cases_by_Neighborhood_joey.csv")) %>% clean_names() 

df_cases=df_c %>%
  filter(as_date(date_reported) == "2022-02-22") %>% 
  separate(neighborhood,into=c("code","name"),sep = ":") %>%
  mutate(code=case_when(
    code=="N35" ~"N0",
    TRUE ~ code
  )) %>%
  select(-objectid,-date_reported)

Regular joining (of dataframes)

neigh2=left_join(neigh,df_cases,by=c("code")) 

tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(neigh2) +tm_polygons("total_positives",alpha=.5)

Joining with other spatial data

Let’s get some data using tidycensus. Need an API key https://api.census.gov/data/key_signup.html



census_api_key("d44395e2fa101f82260fae6b845676d71f017b70")
To install your API key for use in future sessions, run this function with `install = TRUE`.
#what variables
v20 = load_variables(2018,"acs5")
# median_family_income="    B06011_001" 
# all "B00001_001"  
#black "B02009_001"

Get some data:

df_cencus=get_acs(geography = "tract",
                  variables=c("median_inc"="B06011_001",
                              "pop"="B01001_001",
                              "pop_black"="B02009_001"),
                  state="DC",geometry=TRUE,year=2018) 
Getting data from the 2014-2018 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
Using FIPS code '11' for state 'DC'

  |                                                                                
  |                                                                          |   0%
  |                                                                                
  |===============================                                           |  42%
  |                                                                                
  |==========================================                                |  57%
  |                                                                                
  |===============================================                           |  64%
  |                                                                                
  |=================================================                         |  66%
  |                                                                                
  |===========================================================               |  79%
  |                                                                                
  |==========================================================================| 100%
class(df_cencus)
[1] "sf"         "data.frame"
plot(df_cencus)

It’s in long format. Let’s make it wide.

df_cens=df_cencus %>% select(-moe) %>% spread(variable,estimate) 

tm_shape(df_cens) +tm_polygons("median_inc",alpha=.5)

  tm_shape(neigh2) +tm_borders(col="blue",lwd=5,alpha=.2)+
  tm_shape(df_cens) +tm_borders(col="red",lwd=1,alpha=.3)
df_cens_adj=df_cens %>% st_transform(4326)
df_j=st_join(df_cens_adj,neigh2,largest=TRUE)
Warning: attribute variables are assumed to be spatially constant throughout all geometries

Other order?:

Since we want the geometry for the NEIGHBORHOODS, we need a to work a little harder:

df1=df_j %>% select(median_inc,pop,pop_black,objectid) %>%
  group_by(objectid) %>%
  summarise(pop_n=sum(pop),
            pop_black_n=sum(pop_black), 
            adj_median_income=sum(pop*median_inc)/pop_n) 

plot(df1)

#df2=left_join(neigh2,df1)

df2=left_join(neigh2,df1 %>% st_set_geometry(NULL))
Joining, by = "objectid"
df2=df2 %>% mutate(black_perc=pop_black_n/pop_n, covid_rate=total_positives/pop_n)
tm_shape(df2)+tm_polygons(c("adj_median_income","covid_rate","black_perc"))
df2 %>% filter(objectid!=30) %>% tm_shape()+tm_polygons(c("adj_median_income","covid_rate","black_perc"),alpha=.4)

#find where people ride bikes (bikes started in each ‘district’) — Joseph Whites work

bike_data <- read_csv("202209-capitalbikeshare-tripdata.csv") %>% clean_names()
Rows: 386085 Columns: 13── Column specification ────────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): ride_id, rideable_type, start_station_name, end_station_name, member_c...
dbl  (6): start_station_id, end_station_id, start_lat, start_lng, end_lat, end_lng
dttm (2): started_at, ended_at
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
bike_data_sf <- bike_data %>%
  mutate_at(vars(start_lat, start_lng), as.numeric) %>%
  st_as_sf(
    coords = c("start_lat", "start_lng"), 
    agr = "constant",
    crs = "4326"
  ) %>%
  sample_n(1000)
st_crs(bike_data_sf)
Coordinate Reference System: NA
st_crs(neigh2)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
st_crs(bike_data_sf) <- 4326
bike_data_sf_1 <- bike_data_sf %>%
  select(geometry) %>%
  rename(geom_points = geometry)

neigh2_1 <- neigh2 %>%
  select(geometry)
df_plz = st_join(bike_data_sf_1, neigh2_1, join = st_within)

#df_bikes_count <- count(as_tibble(df_plz), )

https://dw-rowlands.github.io/Job_Density_and_Commutes/Job_Density_and_Commutes.html

https://walker-data.com/census-r/index.html

What other interesting questions - is peak in bike data actually people going to work or tourist? - are casual riders in the more touristy areas? open data dc, tourist data, reverse geo coding to find location of monuments :):

Question to answer: Is there a spatial patter driven associated with member/casual? - look on open data dc site - could be some data on tourist activity

LS0tCnRpdGxlOiAiUmVwbGFjZW1lbnQgQ2xhc3MgLSBDbGVhbmluZyB1cCBhbmQgU3BhdGlhbCBKb2lucyIKZGF0ZTogICIyMDIyLTExLTA5IgphdXRob3I6ICJDb2FjaCBTa3VmY2EiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCkFkYXB0ZWQgYW5kIGVkaXRlZCBieSBKb3NlcGggV2hpdGUKCgoKVGhpcyBub3RlYm9va3Mgd2lsbCBoZWxwIG1lIHRhbGsgdGhyb3VnaCBhIGJpdCBvZiBtYXRlcmlhbCByZWdhcmRpbmcgKnNwYXRpYWwgZGF0YSosIAp3aXRoIHNvbWUgaW50cm9kdWN0aW9uIHRvIGBzcGF0aWFsIGpvaW5zYC4KCktleSB0b3BpY3MgZm9yIHRvZGF5OgoKKiBVc2luZyBgc2ZgIHBhY2thZ2UKKiBVc2luZyBgdG1hcGAgcGFja2FnZQoqIFVzaW5nIGB0aWR5Y2Vuc3VzYCBwYWNrYWdlCiogUmVtaW5kZXIgb24gam9pbnMgKGZvY3VzOiBsZWZ0IGpvaW4pCiogT3VyIFNwYXRpYWwgRGF0YQogICAqIE5laWdoYm9yaG9vZHMKICAgKiBKb2luaW5nIHdpdGggbm9uLXNwYXRpYWwgZGF0YQogICAqIENlbnN1cyBkYXRhCiAgICogSm9pbmluZyB3aXRoIHNwYXRpYWwgZGF0YQoqIGlnbm9yZSBodG1sIGluIGdpdCAgKGlmIHRpbWUpCgojIyBQYWNrYWdlcwoKU3RhbmRhcmRzOgoKYGBge3J9CmxpYnJhcnkoa25pdHIpCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGphbml0b3IpCmxpYnJhcnkobHVicmlkYXRlKSAjIGJlY2F1c2Ugd2Ugd2lsbCBwcm9iYWJseSBzZWUgc29tZSBkYXRlcwpsaWJyYXJ5KGhlcmUpICMgYSBwYWNrYWdlIEkgaGF2ZW4ndCB0YXVnaHQgeW91IGFib3V0IGJlZm9yZSB0aGF0IGRvZXNuJ3QgZG8gbXVjaCwgYnV0IC4uLi4KbGlicmFyeShybmF0dXJhbGVhcnRoKQpsaWJyYXJ5KFdESSkKbGlicmFyeSh0aWdyaXMpCmxpYnJhcnkocmdkYWwpCmxpYnJhcnkoc3ApCmBgYAoKU29tZSBhZGRpdGlvbmFsIHBhY2thZ2VzIGZvY3VzZWQgb24gdG9kYXkncyB3b3JrOgoKYGBge3J9CmxpYnJhcnkoc2YpICMgd29ya2luZyB3aXRoIHNpbXBsZSBmZWF0dXJlcyAtIGdlb3NwYXRpYWwKbGlicmFyeSh0bWFwKQpsaWJyYXJ5KHRpZHljZW5zdXMpCgpgYGAKIyMgSW5mb3JtYXRpb25hbCByZXNvdXJjZXMKCiogQW4gb3ZlcmFsbCByZXNvdXJjZSBvbiBtYXBwaW5nIGluIFI6IGh0dHBzOi8vYm9va2Rvd24ub3JnL25pY29oYWhuL21ha2luZ19tYXBzX3dpdGhfcjUvZG9jcy9pbnRyb2R1Y3Rpb24uaHRtbAoqIEEgc3RhcnRpbmcgcG9pbnQgdG8gbGVhcm4gYWJvdXQgYHNmYDogIGh0dHBzOi8vci1zcGF0aWFsLmdpdGh1Yi5pby9zZi9hcnRpY2xlcy8KKiBHZXR0aW5nIHN0YXJ0ZWQgd2l0aCBgdG1hcGA6IGh0dHBzOi8vY3Jhbi5yLXByb2plY3Qub3JnL3dlYi9wYWNrYWdlcy90bWFwL3ZpZ25ldHRlcy90bWFwLWdldHN0YXJ0ZWQuaHRtbAoqIFRoZSBgdGlkeWNlbnN1c2AgcGFja2FnZTogaHR0cHM6Ly93YWxrZXItZGF0YS5jb20vdGlkeWNlbnN1cy9pbmRleC5odG1sCiogVGhlIGJvb2sgb24gYHRpZHljZW5zdXNgIDogaHR0cHM6Ly93YWxrZXItZGF0YS5jb20vY2Vuc3VzLXIvaW5kZXguaHRtbAoKCiMjIFVzaW5nIHRoZSBOZWlnaGJvcmhvb2QgR2Vvc3BhdGlhbCBEYXRhICh1c2luZyAvZGF0YSkKCgpPdXIgZmlyc3QgZGF0YSBzb3VyY2UgY29tZXMgZnJvbSBvcGVuZGF0YS5kYwoKaHR0cHM6Ly9vcGVuZGF0YS5kYy5nb3YvZGF0YXNldHMvRENHSVM6OmRjLWhlYWx0aC1wbGFubmluZy1uZWlnaGJvcmhvb2RzL2Fib3V0CgoKSSB3aWxsIHVzZSB0aGUgR2VvSlNPTiBmaWxlLiAgKE5ld2VyLCBub3QgbmVjZXNzYXJpbHkgYmV0dGVyLCBidXQgLi4uIGEgc2luZ2xlIGZpbGUuICBOb3Qgc21hbGxlciwgYnV0IC4uLiB0aGlzIG9uZSBpcyBub3QgYmlnLikgIAoKCgoKRGF0YSBpcyBlYXNpbHkgcmVhZGFibGUgCmBgYHtyfQpuZWlnaD1zdF9yZWFkKGhlcmUoIkRDX0hlYWx0aF9QbGFubmluZ19OZWlnaGJvcmhvb2RzX2pvZXkuZ2VvanNvbiIpKSAlPiUgY2xlYW5fbmFtZXMoKQpjbGFzcyhuZWlnaCkKYGBgCgpgYGB7cn0KcGxvdChuZWlnaCkKYGBgCgoKCiMjIFJlbWluZGVyIC0gSm9pbnMKCmBgYHtyfQpkZjE9dGliYmxlKGZydWl0PWMoImFwcGxlIiwiYmFuYW5hIiwiY2hlcnJ5IiksY29zdD1jKDEuNSwxLjIsMi4yNSkpCmRmMj10aWJibGUoZnJ1aXQ9YygiYXBwbGUiLCJhcHBsZSIsImNoZXJyeSIsImxlbW9uIiksCiAgICAgICAgICAgZGVzZXJ0PWMoInBpZSIsImNvYmJsZXIiLCJjb2JibGVyIiwiY2hlZXNlY2FrZSIpLAogICAgICAgICAgIGNhbD1jKDQwMCw0MzAsNTAwLDU1MCkpCmRmMQpgYGAKYGBge3J9CmRmMgpgYGAKYGBge3J9CmxlZnRfam9pbihkZjEsZGYyLGJ5PSJmcnVpdCIpCmBgYAoKIyMgSW52ZXN0aWdhdGluZyBqb2luaW5nIHNwYXRpYWwgYW5kIG5vbi1zcGF0aWFsIGRhdGEKCgpDb3ZpZCBjYXNlIGluZm9ybWF0aW9uIGlzIGF2YWlsYWJsZSBmcm9tIG9wZW5kYXRhZGM6CgpodHRwczovL29wZW5kYXRhLmRjLmdvdi9kYXRhc2V0cy9EQ0dJUzo6ZGMtY292aWQtMTktdG90YWwtcG9zaXRpdmUtY2FzZXMtYnktbmVpZ2hib3Job29kL2Fib3V0CgpSZWFkIGNhc2VzIGluZm9ybWF0aW9uOgoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KZGZfYz1yZWFkX2NzdihoZXJlKCJEQ19DT1ZJRC0xOV9Ub3RhbF9Qb3NpdGl2ZV9DYXNlc19ieV9OZWlnaGJvcmhvb2Rfam9leS5jc3YiKSkgJT4lIGNsZWFuX25hbWVzKCkgCgpkZl9jYXNlcz1kZl9jICU+JQogIGZpbHRlcihhc19kYXRlKGRhdGVfcmVwb3J0ZWQpID09ICIyMDIyLTAyLTIyIikgJT4lIAogIHNlcGFyYXRlKG5laWdoYm9yaG9vZCxpbnRvPWMoImNvZGUiLCJuYW1lIiksc2VwID0gIjoiKSAlPiUKICBtdXRhdGUoY29kZT1jYXNlX3doZW4oCiAgICBjb2RlPT0iTjM1IiB+Ik4wIiwKICAgIFRSVUUgfiBjb2RlCiAgKSkgJT4lCiAgc2VsZWN0KC1vYmplY3RpZCwtZGF0ZV9yZXBvcnRlZCkKCgpgYGAKCgojIyBSZWd1bGFyIGpvaW5pbmcgKG9mIGRhdGFmcmFtZXMpCgpgYGB7cn0KbmVpZ2gyPWxlZnRfam9pbihuZWlnaCxkZl9jYXNlcyxieT1jKCJjb2RlIikpIAoKdG1hcF9tb2RlKCJ2aWV3IikKCnRtX3NoYXBlKG5laWdoMikgK3RtX3BvbHlnb25zKCJ0b3RhbF9wb3NpdGl2ZXMiLGFscGhhPS41KQpgYGAKCgojIyBKb2luaW5nIHdpdGggb3RoZXIgc3BhdGlhbCBkYXRhCgpMZXQncyBnZXQgc29tZSBkYXRhIHVzaW5nIGB0aWR5Y2Vuc3VzYC4gIE5lZWQgYW4gQVBJIGtleSAgIGh0dHBzOi8vYXBpLmNlbnN1cy5nb3YvZGF0YS9rZXlfc2lnbnVwLmh0bWwKCgpgYGB7cn0KCgpjZW5zdXNfYXBpX2tleSgiZDQ0Mzk1ZTJmYTEwMWY4MjI2MGZhZTZiODQ1Njc2ZDcxZjAxN2I3MCIpCgojd2hhdCB2YXJpYWJsZXMKdjIwID0gbG9hZF92YXJpYWJsZXMoMjAxOCwiYWNzNSIpCiMgbWVkaWFuX2ZhbWlseV9pbmNvbWU9IglCMDYwMTFfMDAxIiAKIyBhbGwgIkIwMDAwMV8wMDEiCQojYmxhY2sgIkIwMjAwOV8wMDEiCmBgYAoKCkdldCBzb21lIGRhdGE6CgpgYGB7cn0KZGZfY2VuY3VzPWdldF9hY3MoZ2VvZ3JhcGh5ID0gInRyYWN0IiwKICAgICAgICAgICAgICAgICAgdmFyaWFibGVzPWMoIm1lZGlhbl9pbmMiPSJCMDYwMTFfMDAxIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgInBvcCI9IkIwMTAwMV8wMDEiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAicG9wX2JsYWNrIj0iQjAyMDA5XzAwMSIpLAogICAgICAgICAgICAgICAgICBzdGF0ZT0iREMiLGdlb21ldHJ5PVRSVUUseWVhcj0yMDE4KSAKYGBgCgpgYGB7cn0KY2xhc3MoZGZfY2VuY3VzKQpwbG90KGRmX2NlbmN1cykKYGBgCkl0J3MgaW4gbG9uZyBmb3JtYXQuICBMZXQncyBtYWtlIGl0IHdpZGUuCmBgYHtyfQpkZl9jZW5zPWRmX2NlbmN1cyAlPiUgc2VsZWN0KC1tb2UpICU+JSBzcHJlYWQodmFyaWFibGUsZXN0aW1hdGUpIAoKdG1fc2hhcGUoZGZfY2VucykgK3RtX3BvbHlnb25zKCJtZWRpYW5faW5jIixhbHBoYT0uNSkKYGBgCgoKYGBge3J9CgogIHRtX3NoYXBlKG5laWdoMikgK3RtX2JvcmRlcnMoY29sPSJibHVlIixsd2Q9NSxhbHBoYT0uMikrCiAgdG1fc2hhcGUoZGZfY2VucykgK3RtX2JvcmRlcnMoY29sPSJyZWQiLGx3ZD0xLGFscGhhPS4zKQpgYGAKCgoKYGBge3J9CiM8PDw8PDw8IEhFQUQKI2RmX2o9c3Rfam9pbihkZl9jZW5zLG5laWdoMikKIz09PT09PT0KI2RmX2o9c3Rfam9pbihkZl9jZW5zLG5laWdoMixwcmVwYXJlZD1GQUxTRSkKIz4+Pj4+Pj4gYWFmMDFiZTVjZjcyMTgxOWRkMmRmNjE1YWVmN2ExOTk5YmNlYzBjMgpgYGAKCmBgYHtyfQpkZl9jZW5zX2Fkaj1kZl9jZW5zICU+JSBzdF90cmFuc2Zvcm0oNDMyNikKYGBgCgpgYGB7cn0KZGZfaj1zdF9qb2luKGRmX2NlbnNfYWRqLG5laWdoMixsYXJnZXN0PVRSVUUpCmBgYApPdGhlciBvcmRlcj86CgpgYGB7cn0KIzw8PDw8PDwgSEVBRAojI2RmX2pfcmV2ID0gc3Rfam9pbihuZWlnaDIsZGZfY2Vuc19hZGosbGFyZ2VzdD1UUlVFKQojPT09PT09PQojZGZfal9yZXYgPSBzdF9qb2luKG5laWdoMixkZl9jZW5zX2FkaixsYXJnZXN0PVRSVUUpCiM+Pj4+Pj4+IGFhZjAxYmU1Y2Y3MjE4MTlkZDJkZjYxNWFlZjdhMTk5OWJjZWMwYzIKYGBgCgpTaW5jZSB3ZSB3YW50IHRoZSBnZW9tZXRyeSBmb3IgdGhlIE5FSUdIQk9SSE9PRFMsIHdlIG5lZWQgYSB0byB3b3JrIGEgbGl0dGxlIGhhcmRlcjoKCmBgYHtyfQpkZjE9ZGZfaiAlPiUgc2VsZWN0KG1lZGlhbl9pbmMscG9wLHBvcF9ibGFjayxvYmplY3RpZCkgJT4lCiAgZ3JvdXBfYnkob2JqZWN0aWQpICU+JQogIHN1bW1hcmlzZShwb3Bfbj1zdW0ocG9wKSwKICAgICAgICAgICAgcG9wX2JsYWNrX249c3VtKHBvcF9ibGFjayksIAogICAgICAgICAgICBhZGpfbWVkaWFuX2luY29tZT1zdW0ocG9wKm1lZGlhbl9pbmMpL3BvcF9uKSAKCnBsb3QoZGYxKQpgYGAKCmBgYHtyfQojZGYyPWxlZnRfam9pbihuZWlnaDIsZGYxKQoKZGYyPWxlZnRfam9pbihuZWlnaDIsZGYxICU+JSBzdF9zZXRfZ2VvbWV0cnkoTlVMTCkpCgpgYGAKCmBgYHtyfQpkZjI9ZGYyICU+JSBtdXRhdGUoYmxhY2tfcGVyYz1wb3BfYmxhY2tfbi9wb3BfbiwgY292aWRfcmF0ZT10b3RhbF9wb3NpdGl2ZXMvcG9wX24pCnRtX3NoYXBlKGRmMikrdG1fcG9seWdvbnMoYygiYWRqX21lZGlhbl9pbmNvbWUiLCJjb3ZpZF9yYXRlIiwiYmxhY2tfcGVyYyIpKQpgYGAKCgoKYGBge3J9CmRmMiAlPiUgZmlsdGVyKG9iamVjdGlkIT0zMCkgJT4lIHRtX3NoYXBlKCkrdG1fcG9seWdvbnMoYygiYWRqX21lZGlhbl9pbmNvbWUiLCJjb3ZpZF9yYXRlIiwiYmxhY2tfcGVyYyIpLGFscGhhPS40KQpgYGAKCgoKI2ZpbmQgd2hlcmUgcGVvcGxlIHJpZGUgYmlrZXMgKGJpa2VzIHN0YXJ0ZWQgaW4gZWFjaCAnZGlzdHJpY3QnKSAtLS0gSm9zZXBoIFdoaXRlcyB3b3JrCmBgYHtyfQpiaWtlX2RhdGEgPC0gcmVhZF9jc3YoIjIwMjIwOS1jYXBpdGFsYmlrZXNoYXJlLXRyaXBkYXRhLmNzdiIpICU+JSBjbGVhbl9uYW1lcygpCmBgYAoKYGBge3J9CmJpa2VfZGF0YV9zZiA8LSBiaWtlX2RhdGEgJT4lCiAgbXV0YXRlX2F0KHZhcnMoc3RhcnRfbGF0LCBzdGFydF9sbmcpLCBhcy5udW1lcmljKSAlPiUKICBzdF9hc19zZigKICAgIGNvb3JkcyA9IGMoInN0YXJ0X2xhdCIsICJzdGFydF9sbmciKSwgCiAgICBhZ3IgPSAiY29uc3RhbnQiLAogICAgY3JzID0gIjQzMjYiCiAgKSAlPiUKICBzYW1wbGVfbigxMDAwKQpgYGAKCmBgYHtyfQpzdF9jcnMoYmlrZV9kYXRhX3NmKQpzdF9jcnMobmVpZ2gyKQpgYGAKCmBgYHtyfQpzdF9jcnMoYmlrZV9kYXRhX3NmKSA8LSA0MzI2CmBgYAoKYGBge3J9CmJpa2VfZGF0YV9zZl8xIDwtIGJpa2VfZGF0YV9zZiAlPiUKICBzZWxlY3QoZ2VvbWV0cnkpICU+JQogIHJlbmFtZShnZW9tX3BvaW50cyA9IGdlb21ldHJ5KQoKbmVpZ2gyXzEgPC0gbmVpZ2gyICU+JQogIHNlbGVjdChnZW9tZXRyeSkKYGBgCiAgCmBgYHtyfQpkZl9wbHogPSBzdF9qb2luKGJpa2VfZGF0YV9zZl8xLCBuZWlnaDJfMSwgam9pbiA9IHN0X3dpdGhpbikKCiNkZl9iaWtlc19jb3VudCA8LSBjb3VudChhc190aWJibGUoZGZfcGx6KSwgKQpgYGAKCgpodHRwczovL2R3LXJvd2xhbmRzLmdpdGh1Yi5pby9Kb2JfRGVuc2l0eV9hbmRfQ29tbXV0ZXMvSm9iX0RlbnNpdHlfYW5kX0NvbW11dGVzLmh0bWwKCmh0dHBzOi8vd2Fsa2VyLWRhdGEuY29tL2NlbnN1cy1yL2luZGV4Lmh0bWwKCgpXaGF0IG90aGVyIGludGVyZXN0aW5nIHF1ZXN0aW9ucwotIGlzIHBlYWsgaW4gYmlrZSBkYXRhIGFjdHVhbGx5IHBlb3BsZSBnb2luZyB0byB3b3JrIG9yIHRvdXJpc3Q/Ci0gYXJlIGNhc3VhbCByaWRlcnMgaW4gdGhlIG1vcmUgdG91cmlzdHkgYXJlYXM/Cm9wZW4gZGF0YSBkYywgdG91cmlzdCBkYXRhLCByZXZlcnNlIGdlbyBjb2RpbmcgdG8gZmluZCBsb2NhdGlvbiBvZiBtb251bWVudHMgOik6CgpRdWVzdGlvbiB0byBhbnN3ZXI6CklzIHRoZXJlIGEgc3BhdGlhbCBwYXR0ZXIgZHJpdmVuIGFzc29jaWF0ZWQgd2l0aCBtZW1iZXIvY2FzdWFsPwotIGxvb2sgb24gb3BlbiBkYXRhIGRjIHNpdGUKLSBjb3VsZCBiZSBzb21lIGRhdGEgb24gdG91cmlzdCBhY3Rpdml0eQoK